

gen sample = 1 if e(sample) == 1 // Drop all unused observations. This will speed up the margins command significantly
keep if sample == 1

sum advice if sample == 1 // The following generates data for underlying histogram showing # of obs
local min = `r(min)'
local max = `r(max)'
twoway__histogram_gen advice if sample == 1, freq gen(h x) start(`min') width(1)
replace x = . if h == .
sum advice if sample == 1
global obs = `r(N)'
foreach var of varlist h {
	replace `var' = `var' / $obs * 100
	}


estimates restore $model // The following estimates marginal effects at various levels of X
sum advice if e(sample)
local min = `r(min)'
local max = `r(max)'
margins, dydx(treat) at(advice = (`min'(10)`max') age2 = ($xlo $xhi))
matrix margins = r(table)'
svmat margins, names(col)
matrix at = r(at)
svmat at, names(at)
gen n = _n
sum n if at3 != .
local max `r(max)'
forvalues i = 1/`max' {
sum at3 if n == `i'
replace at3 = `r(mean)' if n == `i' + `max' & at3 == .
sum at4 if n == `i'
replace at4 = `r(mean)' if n == `i' + `max' & at4 == .
}
* Prepare stars indicating when conditional marginal effects' p-value is > 0.05
gen stars = b
replace stars = . if pvalue >= 0.05
gen txt = "*" if stars != .
set textsize 100

* Graph
sum at3
local text_x = `r(max)'
foreach i in $xlo $xhi {
sum b if se !=. & at3 == `text_x' & at4 == `i'
local text_`i' = `r(max)'
}

sum h
local yrange2 = `r(max)' * 4

graph twoway bar h x if x < 0, barw(1) lc(red) bc(red) fintensity(100) yaxis(2) ///
		||	bar h x if inrange(x,0,40), barw(1) lc(gold) fintensity(100) bc(gold) yaxis(2) ///
		||	bar h x if x > 40, barw(1) lc(green) bc(green) fintensity(100) yaxis(2) ///
		||	scatter stars at3 if se != . & at4 == $xlo, msymbol(i) mlabel(txt) mlabsize(small) mlabgap(-1.5) mlabposition(11) mlabcolor(black) ///
		||  scatter stars at3 if se != . & at4 == $xhi, msymbol(i) mlabel(txt) mlabsize(small) mlabgap(-1.5) mlabposition(11) mlabcolor(black) ///
        ||  line b at3 if se != . & at4 == $xlo, clwidth(medthick) clpattern(solid) clcolor(black) ///
        ||  line b at3 if se != . & at4 == $xhi, clwidth(medthick) clpattern(solid) clcolor(black) ///
			, xlabel("$xlabel", nogrid labsize(3.5)) ///
            ylabel(, nogrid labsize(3.5) angle(horizontal)) ///
			ylabel(, axis(2) nogrid noticks nolabels) ///
			yscale(noline alt) yscale(range(0 `yrange2') noline alt axis(2)) xscale(noline) legend(off) ///
            yline(0, lcolor(black)) ///
            xtitle("VAA congruence score", size(4)) ///
            ytitle("Treatment effect", size(4) margin(zero)) ///
		    ytitle("" , axis(2) size(4) margin(zero)) ///	
			ysca(axis(2) titlegap(4)) ///
			xsize(6) ysize(4) ///
			title("$title", margin(bottom)) ///
            text(`text_$xlo' `text_x' "$xlow", placement(e) justification(left) margin(left) color(black)) ///
            text(`text_$xhi' `text_x' "$xhigh", placement(e) justification(left) margin(left) color(black)) ///
			scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white) margin(4 $marginright 4 4)) scale(*1.3)
